A new species of lotic breeding salamander (Amphibia, Caudata, Hynobiidae) from Shikoku, Japan

Background Hynobius hirosei is a lotic-breeding salamander endemic to Shikoku Island in western Japan. Significant allozymic and morphological differences have been found among the populations of this species; however, the degree and pattern of intraspecific variation have not been surveyed using a sufficient number of samples. Methods For the taxonomic revision of H. hirosei, we conducted genetic and morphological surveys using samples collected throughout the distribution. Phylogenetic analysis using the cytochrome b region of mitochondrial DNA and population structure analysis using single nucleotide polymorphisms were conducted to evaluate the population structure within the species and the degree of genetic differentiation. Subsequently, a morphological survey based on multivariate and univariate analyses was performed to assess the morphological variation. Results Genetic analyses revealed three genetic groups (Tsurugi, Central, and Nanyo) within H. hirosei, with the Nanyo group distributed allopatrically from the others, and the Tsurugi and Central groups distributed parapatrically with the formation of a hybrid zone between them. The Nanyo group was morphologically distinguishable from the remaining samples, including the topotype of H. hirosei, based on a smaller body size and several ratio values of characters to snout-vent length, longer axilla-groin distance, shorter tail length, shorter internarial distance, longer upper eyelid length, and larger medial tail width. These results support the notion that the Nanyo group is an undescribed species. However, the remaining genetically differentiated groups could not be divided in the present study. Herein, we described the Nanyo group as a new species.


INTRODUCTION
Hynobius Tschudi, 1838 is a salamander distributed in Eastern and Central Asia (Frost, 2022). Hynobius usually resides in water during the larval period, inhabits the forest floor after metamorphosis, and gathers in the water body to breed. Each species has a different habit and belongs to three breeding types: lentic breeding, which involves the deposit of eggs in pools, sulculus, and ponds; lotic breeding, which involves the deposit of eggs in streams; and underground breeding, which involves the deposit of eggs in underground analysis using SNP, and morphological analyses, which revealed fine-scale population structure and intraspecific variation.

Nomenclatural acts
The electronic version of this article in portable document format will represent a published work according to the International Commission on Zoological Nomenclature (ICZN), and hence the new names contained in the electronic version are effectively published under that Code from the electronic edition alone (see Articles 8.5-8.6 of the Code). This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank Life Science Identifiers (LSIDs) can be resolved and the associated information can be viewed through any standard web browser by appending the LSID to the prefix http://zoobank.org/. The LSID for this publication is as follows: urn:lsid:zoobank.org:pub: 1EFFBF1B-1249-4A9F-8AB2-36BFCD6F6B9E. The online version of this work is archived and available from the following digital repositories: PeerJ, PubMed Central and CLOCKSS.

Sample collection
Using the protocol of Nishikawa et al. (2007), we collected salamanders in the field and fully anesthetized them with an acetone-chloroform solution for subsequent processes. Liver and muscle tissues were collected from the anesthetized salamanders and preserved in 99% ethanol or stored in a deep-freezer for molecular analyses. Voucher specimens were then fixed in 10% formalin, and later preserved in 70% ethanol for permanent storage at the Graduate School of Human and Environment Studies, Kyoto University (KUHE). Animal collection and experiments followed the guideline of animal experiments of the university and were approved by the animal experimentation ethics committee at the KUHE (certificate number: 29-A-7, 30-A-7, 20-A-7, 20-A-5).

Genetic sample and DNA extraction
To survey intraspecific genetic variation in H. hirosei, 116 samples covering the entire geographic range of Shikoku were employed (Fig. 1A, Table 1). Total genomic DNA was extracted from the tissues using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany).
The detailed experimental protocols were described by Aoki, Matsui & Nishikawa (2013). Sequences were assembled and checked using the Chromas Pro 1.34 software     H. stejnegeri from Kumamoto Prefecture (KUHE 28011: AB921166) and Salamandrella keyserlingii from Hokkaido Prefecture (KUHE 13057: AB363573) were added. Table 1 lists the citations that contain these sequences. The obtained sequence data were aligned using MAFFT 7 (Katoh & Standley, 2013). Before phylogenetic analysis, we selected the best substitution model using Kakusan4 (Tanabe, 2011) based on the Akaike information criterion (Akaike, 1974) for Likelihood (ML) method and Bayesian information criterion (Schwarz, 1978) for Bayesian inference (BI) method. The ML tree was estimated using RAxML 8.2 software (Stamatakis, 2014). For the confidence of the ML tree, bootstrap analysis was conducted with 1,000 replicates. Nodes with bootstrap values (BS) >70% were considered as sufficiently supported (Huelsenbeck & Hillis, 1993). The BI tree was generated using MrBayes 3.2.6 software (Ronquist et al., 2012). The BI analyses were performed using two parallel runs of four Markov chain Monte Carlo (MCMC) methods for 20 million generations. The first 25% of the generations were discarded as a burn-in, and one of every 100 of the remaining generations was sampled. The convergence of the MCMC runs was checked using TRACER 1.6 software (Rambaut et al., 2014). The uncorrected p-distance among samples was calculated using MEGA 7 software (Kumar, Stecher & Tamura, 2016).

SNP analysis
To estimate the genetic structure of H. hirosei (114 individuals from 26 populations; DDBJ accession number: DRR361514-DRR361627), we used genome-wide SNP obtained based on MIG-seq (Suyama & Matsuki, 2015) using the Illumina MiSeq system. We followed the protocol described by Suyama & Matsuki (2015), except that an annealing temperature of 38 C (originally, 48 C) for the first PCR and indexed forward and reverse primers for the second PCR (originally, common forward and indexed reverse primers) were employed. The sequenced data were filtered using fastp ver. 3 (Chen et al., 2018) to trim the primer regions and remove low-quality and short reads. The qualified quality phred, length required, and front/tail trimming settings were 30, 80, and 14, respectively. After filtered reads 1 and 2 were connected to one, Stacks ver. 2.3c (Catchen et al., 2011) was executed for SNP detection using the 'denovo_map.pl' program as follows: number of mismatches allowed between stacks within individuals (M) = 2; and number of mismatches allowed between stacks between individuals (n) = 2. Of the several stages in the program, 'ustacks' was configured as follows: minimum depth of coverage required to create a stack (m) = 3 and maximum distance allowed to align secondary reads to primary stacks (N) = 2. Finally, the 'populations' program was set as follows: minimum percentage of samples in a population (r) = 0.7; minimum number of populations in a locus (p) = 2; and specify a minimum minor allele frequency required to process a SNP (min-maf) = 0.05; and specify a maximum observed heterozygosity required to process a SNP (max-obs-het) = 0.75.
The genetic structures within H. hirosei were estimated using Bayesian clustering in STRUCTURE ver. 2. 3. 4 (Pritchard, Stephens & Donnelly, 2000), as well as EasyParallel (Zhao et al., 2020), which is multithreaded parallelization to facilitate population structure analysis. The Monte Carlo Markov chains ran for 1 million generations, including a burn-in of 100,000 generations. This calculation was repeated ten times for the number of each cluster to evaluate genetic differentiation in nuclear markers among the clusters. Multiple STRUCTURE runs were combined using CLUMPP (Jakobsson & Rosenberg, 2007).
When the groups recognized in H. hirosei based on molecular analysis were compared (see Results), only adult males were used as H. hirosei is known to exhibit sexual dimorphism in most of morphological characters (Nishikawa et al., 2007) and a larger number of adult males than adult females was obtained. Sex and maturity were determined by direct gonad observations. To examine differentiation between the groups, we compared SVL using the Student's t-test. The Mann-Whitney's U-test was used to assess meristic characters and ratio values of metric characters to SVL. To examine the overall morphological variation between the groups, linear discriminate analysis (LDA) was conducted using log e -transformed metric values of the total of the 26 measurements. The difference in the LDA score between recognized groups was then determined using the Mann-Whitney's U-test.
The clutch size among the populations was compared using Tukey's pairwise post-hoc test.
In these analyses, the TAL data for the regenerated tail were omitted. All statistical analyses were performed using Past 4.04 (Hammer, Harper & Ryan, 2001).
For larval specimens, the developmental stage was determined according to Iwasawa & Yamashita (1991) and the SVL was measured accordingly.

Mitochondrial DNA analysis
Complete sequences of the cyt b (1,141 bp) of mtDNA were obtained, which included 450 variables and 345 parsimonious informative sites. The best selected substitution models were the codon-equal rate model with the general time reversible model (GTR: Tavaré, 1986) + G for the first, second, and third in the ML analysis and the codon proportional model with SYM_Gamma (Zharkikh, 1994), HKY85_Gamma (Hasegawa, Kishino & Yano, 1985), and GTR_Gamma for the first, second, and third in the BI analysis, respectively. MCMC analysis using MrBayes did not converge; thus, only the phylogenetic tree constructed by ML (likelihood value [lnL] = −7,206.929; Fig. 1B The Nanyo group had a sister relationship with H. sematonotos (BS = 99). Individuals of both the Tsurugi and Central groups were sympatrically found at locality 10.
The mean pairwise genetic distances (Table 2) were 9.4% between the Tsurugi and Central group, 10.3% between the Tsurugi and Nanyo group, and 6.9% between the Central and Nanyo group. The genetic distance between the Nanyo group and H. samatonotos was 4.6%.

SNP analysis
A total of 373,642 reads generated based on MIG-seq were filtered using fastp, resulting in 287,094 reads. Stacks detected 716 SNP loci in the reads, and the genotyping rate was 0.74. Based on STRUCTURE analysis, at K = 2, individuals were clustered into southeastern Shikoku (blue) vs. the remaining area of Shikoku (orange), and transitional structures were found at localities 8 and 9 (Sanuki Mountains) and 10-12 (Fig. 2). Locality 7 was included in the Central group by mitochondrial analysis, but was judged as the Tsurugi group in the STRUCTURE analysis. At K = 3, individuals from southeastern Shikoku (blue), the remaining area of Shikoku, except the Onigajo and Sasayama Mountains (orange), and the Onigajo and Sasayama Mountains (localities 24-26; green) were divided. These three groups recognized by mtDNA analysis were also detected using the SNP data; however, the Tsurugi and Central groups were not separated and clinally continued in the SNP data. Thus, we temporarily treated these two groups as one group for the morphological analyses.

Morphological analysis
To survey the morphological separation of the Nanyo group from the other groups, we conducted morphological analyses of the Nanyo and Tsurugi + Central groups.
The obtained values of SVL and the ratios of the other characters to the SVL are shown in Table 3   In the LDA results, the eigenvalues of the axis accounted for 2.79. On the axis, the highest absolute magnitude of the standardized discriminant coefficients was −266.35 of SVL, followed by that of TRL (224.66), HL (72.85), AGD (−27.06), and BTAW (−10.84). The Nanyo group and the remaining samples were significantly separated (U = 0.00, P = 0.0000), with a small overlap observed (Fig. 3).

Systematics
In summary, the Nanyo group and the remaining samples were clearly separated based on the nuclear genome and external morphology, whereas the Tsusrugi and Central groups formed hybrid zones. The Nanyo group and the remaining samples were genetically isolated; thus, the Nanyo group must be an independent evolutionary lineage with sufficient genetic and morphological differentiation. This information strongly indicates that the Nanyo group is an independent species. The remaining samples (hereafter H. hirosei sensu stricto) included topotypic H. hirosei from Mt. Ishizuchi; however, the Nanyo group has not been named. Therefore, the Nanyo group was revealed to be a new species.
Diagnosis: A large-sized species (adult SVL 73.6-87.5 mm in males) of the lotic-breeding Hynobius, breeding in montane streams; dorsum uniformly dark reddish brown and immaculate in adult; tips of fore-and hindlimbs adpressed on body scarcely meeting (overlap of −2.0 to 0.0 costal folds in males); fifth toe well developed; ova large, pigmentless; egg sacs relatively long and crescent in shape, with distinct whiptail structure on free end; larvae lack claws on their tips of fingers and toes; most similar to H. hirosei, but distinct based on its smaller body size, longer axilla-groin distance, shorter tail length, shorter internarial distance, longer upper eyelid length, and larger medial tail width. Hynobius oni is genetically closer to H. sematonotos than H. hirosei based on mtDNA; however, H. oni has no large markings on the body, in contrast to many silvery spots on H. sematonotos, and a larger SVL than H. sematonotos.
Description of holotype: Head-body moderately large and robust; head oval and moderately depressed, distinctly longer than wide; snout rounded, slightly projecting beyond lower jaw; nostril close to snout tip; labial fold absent; eye large, prominently protruded, slightly inset from edge of head in dorsal view; upper eyelid well developed, shorter than snout; gular fold distinct, curving slightly anteriorly; parotoid gland evident, extending from angle of jaw to gular fold; postorbital grooves distinct, branching posterior to angle of jaw, one short and running down to lower jaw, the other long and posteriorly to parotoid gland; vomerine tooth series wider than long, V-shaped, anterior margin distal to choanae; tongue broad, both sides free from mouth floor; fore-and hindlimbs long and thick; number of costal grooves between axilla and groin 13; depressed limbs separated by two costal folds; relative length of fingers I<IV<III<II, toes I<V<II<IV<III; fifth toe well developed; cloaca longitudinal slit; genital tubercle on anterior cloaca absent; tail moderately short and thick, cylindrical at base, increasingly compressed posteriorly, dorsal fin evident posteriorly; tip of tail rounded in lateral view. Color: Dorsum uniformly dark reddish brown without marking; underside lighter than the dorsum; underside of tail slightly ochre; iris dark brown without marking. In preservatives, dorsal coloration tends to fade and becomes gray-brown, but otherwise has no obvious changes.
Variation: Morphometric data are presented in Table 3. Individuals of the type series and referred specimens are generally similar in body size and proportion.  ; df = 110, F = 16.01, P = 0.0000) and from the Tsurugi Mountains, Tokushima Prefecture (mean: 31.0 ± 12.4, range: 15-52, n = 9 (Tamura, 2012); df = 110, F = 16.01, P = 0.0356). Egg sacs are crescent in shape (length 176.8 ± 23.8 mm (n = 15), width 17.0 ± 1.9 mm (n = 15)) with relatively thicker envelope than that of other congeners, but thinner than that of H. boulengeri, and similar in thickness to that of H. naevius, H. oyamai, and H. sematonotos, with a distinct whiptail structure on the free end (length 35.2 ± 23.7 mm (n = 15)). Both the animal and vegetal poles of eggs have a cream color. The eggs (oval diameter: 4.9-5.9 mm) form a single and/or double row in each egg sac.
Larvae: SVL for fully grown larvae at stage 63 (n = 3) was 21.7 ± 1.5 mm, and SVL for a larva when onset of metamorphosis in late July at stage 64 (n = 1) was 20.6 mm; head rounded in dorsal and lateral views (Fig. 6); snout short and broadly rounded; eyes slightly protruded, inset from the edge of head in dorsal view; labial fold distinct at posterior half of upper jaw; external gills developed; caudal fin higher than head; dorsal fin higher than ventral fin; origin of dorsal fin at distal half to three-fourth of trunk; ventral fin originating from vent; tail tip weakly pointed; limbs slender; claws on fingers and toes absent. In life, dorsum light brown with small dark-brown dots and blotches; venter whitish and transparent; golden dots scattered on tail fin. In preservative, the dorsal coloration tends to fade, becoming light brown and golden dots fading to white.  1A).
Natural history: Hynobius oni inhabits around mountain streams with partially exposed bedrock (Fig. 8). On April 29, 2019, egg sacs, some females with egg, and many males were observed in the water. On April 24, 2021, 12 adult males and egg sacs were found under stones in the water. Thus, the breeding season is presumed to be late April. Adults in the non-breeding season and juveniles were found under stones, rotten wood, and debris near the stream. However, the larval life history of H. oni is poorly understood. Only one overwintered larva (total length: 71 mm) was once found in an open stream through three-night surveys on April 29, May 3, 2019 and April 28, 2022. In the surveys, the observed density of overwintered larvae was lower (approximately 0.25 per 1 m 2 ) than that of H. hirosei (more than 10 per 1 m 2 ; our personal observation) on May 8, 2021. No sympatric hynobiid salamanders (Hynobius or Onychodactylus) were observed with H. oni.
Conservation: As the known range of the new species is restricted to small areas of only the Onigajo and Sasayama Mountains, habitat destruction and over-collecting for private purposes could have negative impacts on the natural populations, as reported for Onychodactylus tsukubaensis (Yoshikawa & Sakamoto, 2016), which also occurs in small areas (the Tsukuba Mountains of Ibaraki Prefecture). Hynobius oni has a very small range and some habitats have been degraded; thus, this species should be protected as an endangered species. Furthermore, the conservation status of H. hirosei should be reconsidered immediately after this taxonomic change (presently, listed as Near Threatened (Matsui, 2014)).

DISCUSSION
Among the three groups of H. hirosei sensu lato observed in this study, the Nanyo group (Hynobius oni) was demonstrated to be an independent species based on genetic and morphological evidence. However, regarding the Tsurugi and Central groups, whether they are distinct species remains unclear until the nature of the hybrid zone is revealed, as they could not be genetically isolated from each other in the SNP analysis. Further analyses are necessary to determine the taxonomic status of the two groups by estimating the degree and demography of the hybrid zone, and elucidate the process of separation of the groups using robust DNA data, as reported by Burbrink & Ruane (2021).
Terrestrial animals inhabiting Shikoku tend to show genetically high inter-and intraspecific divergence (Kato, Morii & Tojo, 2013;Dejima & Sota, 2017;Tominaga et al., 2019b). The three groups with high divergence in H. hirosei sensu lato were distributed in different mountainous areas, and their distributions were usually separated by large rivers and lowlands. On the other hand, hybrid zones between the Tsurugi and Central groups were found, where mountains continued without being divided by large rivers. Large rivers are suggested to contribute to the distribution and maintenance of the population in H. hirosei sensu lato. Furthermore, the same pattern of distribution as the three groups in H. hirosei sensu lato was found in the harvestman, Pseudobiantes japonicus (Kumekawa et al., 2014(Kumekawa et al., , 2019, which is also a terrestrial and low-dispersal animal, similar to the salamander analyzed in this study. Therefore, the concordant pattern recognized in the salamander and harvestman suggests that large rivers may cause geographic variations in Shikoku. Hynobius oni has a significantly smaller body size than H. hirosei sensu stricto. Based on our preliminary survey, H. oni tended to have a smaller number of overwintered larvae and possibly a smaller metamorphosing size than H. hirosei. In the Onigajo and Sasayama Mountains, where H. oni occurs, smaller numbers of stable streams with deep pools were found relative to other areas in Shikoku as the steep riverbeds were formed by exposed bedrock. Some streams dried up during winter. In such shallow and unstable water levels, larvae of lotic-breeding Hynobius do not overwinter often and do not grow to a large metamorphosing size (H. kimurae in Misawa & Matsui (1997); H. osumiensis in ), which agree with our preliminary field surveys in the Onigajo and Sasayama Mountains. Additionally, as shown in other congeners (Misawa & Matsui, 1997), a smaller metamorphosing size may induce smaller adult size in H. oni.
Terrestrial salamanders are well known to coexist through a partition of food resources based on their body size difference (Vignoli, Bissattini & Luiselli, 2016). In the Onigajo and Sasayama Mountains, only H. oni occurs. The smaller adult size of H. oni (79.6 mm in mean SVL) relative to that of H. hirosei sensu stricto could not permit its co-occurrence with other syntopic smaller-sized salamanders from Shikoku (e.g., H. kuishiensis with 60.8 mm, H. tsurugiensis with 62.0 mm (Tominaga et al., 2019b), and Onychodactylus kinneburi with 72.6 mm in mean SVL (Yoshikawa et al., 2013)).
In contrast, H. hirosei sensu stricto is a large species of the genus. Nishikawa et al. (2007) first reported that species of the H. boulengeri complex (including H. hirosei sensu stricto) tended to have a larger body size when they occurred sympatrically with other hynobiid salamander. Such interspecific morphological differentiation is known to occur to avoid interspecific competition for resources, such as food and habitat (Brown & Wilson, 1956;Adams & Rohlf, 2000;Melville, 2002). Indeed, H. hirosei sensu stricto occurs syntopically with H. kuishiensis, H. tsurugiensis, and O. kinneburi. Thus, this sympatric distribution might be possible due to the body size differences between them (Nishikawa et al., 2007).
The range of H. oni is one of the smallest among lotic breeding salamanders in Shikoku. Currently, this range has reduced owing to loss of habitat. Loss and fragmentation of habitats have the most negative impact on amphibian populations and their genetic diversity (Reh & Seitz, 1990;Frankman, Ballou & Briscoe, 2002;Cushman, 2006), especially species that have small distributions, such as H. oni. Presently, local populations of H. hirosei sensu lato are designated as CR+EN in the Red Data Book (RDB) of Matsuyama City (Okayama, 2012) and EN in the RDB of Kagawa Prefecture (Shinohara, 2021). Further, habitats in Kagawa Prefecture are limited and their population sizes are extremely small compared to those in other areas. The recent pet trade of rare salamanders have become a serious conservation problem (Terui & Tokuda, 2021). Therefore, legal protection of H. oni and H. hirosei sensu stricto by national and local governments is urgently needed.

CONCLUSIONS
In the present study, the actual genetic population structure and degree of genetic divergence within Hynobius hirosei, which has been reported to have large genetic intraspecies divergence, were evaluated using mtDNA and nuDNA markers (SNP). Phylogenetic analysis using mtDNA revealed three divergent lineages, including the Tsurugi, Central, and Nanyo groups (genetic distance: 6.9-10.3% in cyt b). Further, STRUCTURE analysis using SNP revealed that the Nanyo group is genetically isolated from the other groups and the Tsurugi and Central groups form hybrid zones. Morphological analyses also revealed that the Nanyo is distinct from the other groups. Collectively, these results strongly indicate that Nanyo group is a distinct species, and is referred to as H. oni sp. nov. in this study.